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We study the dynamical response of triangular Josephson junction arrays, modelled as a network 
of resistively- and capacitively-shunted junctions (RCSJ's). A flux flow regime is found to extend 
between a lower vortex-depinning current and a higher critical current, in agreement with previous 
calculations for square arrays. The upper current corresponds either to row-switching events accom- 
panied by steplike jumps in the array resistance, or to a depinning of the entire array. In the flux 
flow regime, the dynamical response to the bias current is roughly Ohmic, and the time-dependent 
voltage can be well understood in terms of vortex degrees of freedom. The vortex friction coefficient 
rj depends strongly on the McCumber-Stewart parameter /3, and at large /3 is approximately inde- 
pendent of the shunt resistance R. To account for this, we generalize a model of GeigenmuUer et al 
to treat energy loss from moving vortices to the phase analog of optical spin waves in a triangular 
lattice. The value of rj at all values of (3 agrees quite well with this model in the low-density limit. 
The vortex depinning current is estimated as 0.0427c, independent of the direction of applied cur- 
rent, in agreement with static calculations by Lobb et al. A simple argument suggests that quantum 
effects in vortex motion may become important when the flux flow resistivity is of order h/{2e)^ per 
unit frustration. 

PACS numbers: 74.50. -)-r, 74.60. Ge, 74.60.Jg, 74.70.Mq 
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I. INTRODUCTION 

The behavior of vortices in Josephson junction arrays 
(JJA's) has attracted much recent attention. Such 
vortices are coherent arrangements of phases of the super- 
conducting order parameter, which may move through 
the array hke particles, in response to forces generated 
by external currents. They can be generated by a mag- 
netic field, or excited thermally. At low velocities, it has 
been proposed that the motion of a vortex in a square 
array can be described by a Josephson-like equation of 
the form^ 

(1) 

Here x is the vortex position along a line through the 
plaqucttc centers and perpendicular to the external cur- 
rent /, a is the lattice constant, R and C arc the shunt 
resistance and shunt capacitance of each junction, and Id 
is the vortex depinning current, which is linearly related 
to the junction critical current Ic- This equation may 
describe the behavior of a vortex in a square array at low 
velocities. 

Because of the mass term in this equation, a vortex 
might be expected to move ballistically under appropri- 
ate circumstances - that is, a vortex, once set into mo- 
tion, would remain in motion even if the driving current 
is turned off. There are several experimental reports of 
such motion*'^. However, numerical calculations and an- 
alytical studies, based on classical equations of motion, 
have not yet produced such ballistic motion^''^'^^'^"'. The 
numerical studies have been carried out either in square 
arrays^'^'^^, or using a simplified representation of the 
Josephson interaction, in triangular arrays''. Instead of 
ballistic motion, in these simulations, it is generally found 
that the junctions in the wake of the vortex oscillate, 
causing the vortex to lose its energy to the array. In an 
analytical study based on a continuum model^°, an equa- 
tion of motion for individual vortices is derived for both 
square and triangular arrays, starting from an effective 
action for the array. This paper concludes that a narrow 
window of vortex velocity exists in a triangular array, 
for which ballistic motion may be possible. However, it 
seems unlikely that this window is the regime which is 
probed experimentally^'® 

In this paper, we present dynamical calculations for tri- 
angular Josephson junction arrays. The calculations are 



carried out within a classical model of rcsistivcly- and 
capacitively-shunted junctions (RCSJ's). Our principal 
aim is to study the motion of single vortex in the pres- 
ence of an applied d. c. current, to see if ballistic vortex 
motion is possible. Since the depth of the vortex poten- 
tial in a triangular array is believed to be smaller than 
that of a square array such motion would seem more 
likely, at first glance, than in square arrays. However, 
we find no such ballistic motion under any circumstances 
investigated. Instead, the motion of the vortex, when it 
is a coherent excitation of the array, generally falls into a 
"flux flow" regime, where the vortex moves with approx- 
imately Ohmic resistivity, described by a characteristic 
vortex viscosity. 

As in previous simulations in square arrays'''^'^"^'^'*, and 
in similar studies in triangular arrays based on a sim- 
plified Josephson coupling^, we find that the flux flow 
region terminates at sufficiently high current in either 
of two ways. One possibility is for the entire lattice of 
Josephson junctions to be "depinned" causing the volt- 
age to increase sharply. A second decay mode, which 
predominates in sufficiently underdamped arrays, is for 
the vortex to excite row switching events, which pro- 
duce steplike increases of the array resistance as an entire 
row of junctions switches from superconducting to nor- 
mal. Both such decay modes have been in a variety of 
experiments^'^' 

Perhaps the most striking result of our simulations, 
also reported previously in square arrays^'^^, is the per- 
sistence of the quasi-ohmic flux-flow regime even in high- 
resistance arrays where, according to simple models, bal- 
listic motion should be possible. To account for this, 
we generalize a model of Geigenmiiller et aP to obtain 
a simple, nearly analytic model from which an effective 
vortex friction coefficient can be computed at any value 
of the vortex velocity and junction McCumber-Stewart 
coefficient^® 

P = 2eB?IcClh. (2) 

(3 is a dimcnsionless measure of damping in a single junc- 
tion (large j3 means low damping). In our generaliza- 
tion, we take explicit account of the array lattice struc- 
ture, so that a short-wavelength cutoff appears naturally 
in the resulting expression for the vortex friction coeffi- 
cient. The model gives semiquantitative agreement with 
our numerical experiments. In particular, like the contin- 
uum model of Ref. 9 for square arrays, it accounts for the 
persistence of this friction coefficient even at high values 
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of the shunt resistance R. The present model is, in prin- 
ciple, applicable to dissipation from an arbitrary density 
of vortices moving with arbitrary velocities, and to losses 
produced by a. c. external currents. However, we test it 
here only for single vortex motion. 

The remainder of this paper organized as follows. In 
Section II, we describe our calculational model and nu- 
merical method. Section III reports the results of our 
critical current calculations in triangular arrays both 
with and without vortices. Section IV describes in 
greater detail our numerical results in the flux flow 
regime, and Sections V and VI summarize the vortex 
friction model, and our investigation of possible ballistic 
vortex motion. A brief discussion follows in Section VII. 

II. MODEL 

The relevant geometry of our triangular array is shown 
in Fig. 1. Within the array interior of the array, each su- 
perconducting grain is connected to its six nearest neigh- 
bors via Josephson coupling. The boundary conditions 
involve fixed external current injection. We consider two 
directions of current injection, as illustrated in the Fig. 
1, the so-called [101] direction, and the [211] direction^". 

We describe the dynamics within the RCSJ model 
at zero temperature, as previously described for square 
arrays^^'^^. In this model, the current through each junc- 
tion is the sum of three terms: a charge flow through 
an effective intergrain capacitance, a current through a 
shunt resistance, and a Josephson supercurrent. We as- 
sume that the supercurrent is sufficiently weak that we 
can discard the induced screening current. 

With these assumptions, the current between grains i 
and j is: 

Here lij is the total current from grain i to grain j; Cij 
and Rij are the shunt capacitance and shunt resistance 
between grain i and grain j; I^-ij is the critical current of 
the Josephson junction between grain i and grain j; (pi 
is the phase of the order parameter on grain i. Vij and 
Aij are the voltage difference and magnetic gauge phase 
factor between grain i and grain j, defined by 

Vij^Vi-Vj = ^^f^{cl^i-^,), (4) 

and 



where $0 = he/ (2e) is the flux quantum, A is the vec- 
tor potential of the applied magnetic field and Xi is the 
position of the center of grain i. In the present paper, 
we include only the intergrain capacitance, discarding the 
capacitance between the grains and ground. Current con- 
servation at each grain is described by Kirchhoff 's Law, 

i 

where li-ext is the external current fed into grain i. We 
assume that all the capacitances, critical currents, and 
shunt resistances have unique values C, Ic, and R. Fi- 
nally, the use of classical equations of motion implies the 
assumption that quantum effects^^ arising from the non- 
commutativity of charge and phase variables can be ne- 
glected. 

Our boundary conditions are shown in Fig. 1. In the 
direction of current injection, we introduce a uniform 
current li-^ext = I into each boundary grain along one 
edge, and extract the uniform current from each bound- 
ary grain on the other edge. In the transverse direction, 
we use periodic boundary conditions, as in our previous 
work^**. The gauge factor Aij satisfies 

plaquette 

where B is the magnetic field strength; S is the area of 
each triangular plaquette; $0 is the flux quantum. / 
is the so-called frustration. Note that the primitive cell 
of a triangular array, consists of two adjacent triangular 
plaquettes. 

The coupled equations (3), (4) and (5) are solved as 
described previously We use a fourth-order Runge- 
Kutta algorithm with time step At, where At ranges 
from O.Olto to 0.05to (^o = h/{2eRIc) is a characteris- 
tic damping time) , depending on the desired precision of 
calculation. Further details may be found in Ref. 14. 

III. CRITICAL CURRENT AND VORTEX 
DEPINNING CURRENT 

When / = 0, with d. c. bias current injected in the 
[101] direction, we find a critical current of exactly 21^. 
This value can be easily understood, since in this case no 
current passes through the junctions perpendicular to the 
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bias current, so that the entire array behaves much hke a 
single junction. With [2ll] current injection, the critical 
current is found to be approximately 1.76/c. This value 
can be understood by considering the single triangular 
plaqucttc shown in Fig. 2. For such an arrangement, 
the injected current / is related to the phase difference 
4> by I / Ic = s'm(j> + sin(2(^). The right hand side cannot 
exceed {I/Ic)max = 1.7602, which corresponds to our 
numerically obtained array value. We have checked the 
phase configuration for each grain in the array in this 
geometry, and find that it decomposes exactly into unit 
cells of this type. 

Next, we discuss the dynamical response of an array 
under d. c. bias and in the presence of a single vortex. 
Such a vortex can be introduced by considering an N x N 
array containing 2^^^ triangular plaquettes, and a flux 
/ = 1/{2N'^) {N being the number of junctions spanning 
the array in one direction, as in Fig. 1). Table 1 lists the 
critical currents Id for / = 1/{2N'^) as a function of N for 
[101] current injection. This critical current appears to 
be independent of /3, at least in the range < /3 < 1000. 
The critical currents are extracted from an J — (V) plot, 
such as shown in Fig. 3 for [101] current injection aX (3 — Q 
and (3 = 10, and Fig. 4 for [211] direction at /? = 0. 
Extrapolating by eye a plot of Id versus 1/A'' towards 
N = oo at f = 1/{2N'^), we estimate a critical current 
of about 0.042/c for bias current injected in the [101] di- 
rection. In the [2ii] direction, for = 8, we estimate 
a value of about 0.0417c, possibly dependent on the ini- 
tial phase conflguration. In both cases, our calculated 
critical currents are in reasonable agreement with those 
obtained by Bobbert'^, using a piece- wise linear function 
to approximate the sinusoidal coupling function. 

The array critical current at fleld / = 1/{2N'^) can be 
interpreted as the depinning current of a single vortex. 
It is of interest to compare our calculated value with the 
energy barrier for depinning. This energy barrier was 
calculated by Lobb et al,^^, who used static methods to 
obtain a value of about 0.043?i7c/(2e). This represents 
the energy which must be overcome in order to move a 
vortex from the center of one triangular plaquette to the 
center of an adjacent plaquette. 

To make this comparison, and also to account for 
the apparent isotropy of the vortex depinning current, 
we have used a simple model for the vortex potential 
U{r)[r = {x,y)]. Since U{x,y) must have the array pe- 
riodicity, we express it as a Fourier series involving only 
Fourier components from the reciprocal lattice. The sim- 



plest approximation consistent with the point symmetry 
of the lattice is to include only the smallest-magnitude 
Fourier components, i. e. 

f/(r) = Uo + UiY^ cos(K • r) (6) 

K 

where we take Ui > and K is one of the six smallest 
reciprocal lattice vectors, (r = is interpreted as a grain 
center, and hence a maximum in the potential.) The 
potential barrier for vortex depinning in this model is 
readily shown to be just Ui. 

To estimate the critical current in this picture, we add 
to the vortex potential energy a term ]$oJ x r|/c, where 
J is the external current density. This term corresponds 
to the Magnus force J x z$o/c on a single vortex. When 
this term is included, we find numerically that the barrier 
for vortex motion between adjacent triangular plaquettes 
disappears at ahJ/{2e) k, 0.984J7i for current injected 
in the [101] direction. Taking Ui = QM'iMc/ {2e) from 
the results of Lobb et af^, and using J = I/a for [101] 
current injection, we see that our calculated vortex de- 
pinning current of 0.042Jc is in excellent agreement with 
the static results. This calculation also suggests that ex- 
pression (6) is a reasonable approximation for the vortex 
potential. 

IV. SINGLE VORTEX MOTION: NUMERICAL 
RESULTS 

As shown in Fig. 3 and Fig. 4, the / — (V) character- 
istics display a long low voltage tail at currents above 
the vortex depinning current. This current regime is ap- 
proximately Ohmic. Since the / — {V) characteristics in 
this region can be understood in terms of single vortex 
motion in the array, this region is often called the flux 
flow regimc^^^. 

At /3 = 0, the flux flow regime extends to about 1.9/c 
and 1.5/c for current biased in the [lOl] and [2ll] direc- 
tions. Both of these values are close to the critical current 
of the array at / = 0. Above these currents, the whole 
array is depinned and the vortex picture is not applica- 
ble. At sufficiently large values of (3, with the current 
injected at [101] direction, the flux flow regime termi- 
nates at lower currents, where there is a "row switching" 
event^^'"'^^ rather than the depinning of the entire array ~ 
that is, one or more rows of junctions parallel to the di- 
rection of current injection switch from the supercurrent 
state to the resistively dissipative state. This occurs, for 
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example, near / = 1.5/c at /? = 10 in 8 x 8 array. Above 
this row-switching threshold, the picture of ohmic resis- 
tance by flux flow of vortices is no longer valid. At yet 
higher /3 values, the flux flow regime terminates at smaller 
currents, and there may be more than one row switching 
event before the entire array is depinned. At (3 — 100 
in 8 X 8 array, for example, we find two row switching 
events in our calculations. As in Ref. 9, we also find the 
staircase-like structure in the / — {V) curve in the flux 
flow region, which may arise from the interaction of the 
vortex with its image neighboring vortices generated by 
the periodic boundary conditions. 

Figs. 5(a) and 5(b) show the time-dependent space- 
averaged voltage drop V{t) across the array - that is, 
the difference between the average voltage on the line of 
grains where the current is injected, and the line from 
which it is extracted - at two current values in the flux- 
flow regime. In both cases, {3 = and the current is 
injected in the [101] direction. V{t) is characterized by 
periodic sharp peaks which resemble the time-dependent 
voltage of a single junction, the frequency of which in- 
creases with increasing bias current. 

If the single vortex picture is correct, the spike fre- 
quency is related to (V). The period of oscillation 
should correspond to the motion of a vortex by one unit 
cell. Furthermore, a complete vortex circuit around the 
lattice should produce a phase change of 2tt across the 
array. Then using the Josephson relation, we obtain 

where T is the period for one complete vortex circuit. For 

an N X N array, T = N/vy. Our numerical results are 
in excellent agreement with this relation, thus confirming 
the vortex motion picture in the fiux fiow regime. 

By examining the time-dependent voltage of each sin- 
gle junction, one can also deduce the actual vortex path 
in the array. This path is displayed in Fig. 5(c) for an 
8x8 array with bias current in the [lOl] direction. We 
find that this path is independent of current magnitude 
in the [101] direction. As shown in the figure, it is a 
straight trajectory through the middle of the array. 

The behavior of vortices is more complex when current 
is applied in the [211] direction. In this case, the flux flow 
regime in an 8 x 8 overdamped array consists of three dis- 
tinct subregions with different slopes as shown in Fig. 4. 
A typical voltage trace V{t) from each subregion is shown 
in Figs. 6(a), 6(b), and 6(c). The relationship between 
spike frequency and time-averaged voltage still holds, im- 



plying that the picture of vortex motion is still correct 
in this direction. However, the vortex path is different 
in each of the three subregions. These paths, as deduced 
from the time-dependent voltages of each junction, are 
shown in Figs. 6(d), 6(e), and 6(f). 

At sufficiently low bias current in [101] direction, V{t) 
shows a double-peaked structure. We believe that this 
structure originates in the special geometry of the tri- 
angular lattice, in which each primitive cell has two tri- 
angular plaquettes which are inequivalent. When a bias 
current is applied in the [101] direction, the vortex will 
pass alternately through each of these plaquettes, some- 
how producing a double-peaked structure in V{t). As 
the bias current increases, the double-peaked structure in 
V{t) seems to disappear. However, the time-dependent 
voltage of each individual junction still exhibits a double- 
peaked structure. The disappearance of the double-peak 
structure in V{t) is therefore due to space-averaging. We 
conclude that our simple vortex potential is qualitatively 
correct even at higher bias current. Of course, above the 
array depinning current (/ = 27c for this direction), the 
vortex picture breaks down and V{t) shows no simple be- 
havior, just as was found previously for square arrays^'^^. 

V. FLUX FLOW RESISTIVITY AND VORTEX 
FRICTION COEFFICIENT 

As noted above, in the flux flow regime, the Joseph- 
son network behaves approximately ohmically. In this 
ohmic regime, we can define a vortex friction coefficient 
7] by equating the driving force J$o/c to the frictional 
drag force riv, where v is the vortex velocity. This gives 
(assuming current injected in the [101] direction) 



where / — aJ is the current injected into one node. 

This vortex friction coefficient can be estimated in a 
simple way by c;quating the frictional losses to the power 
dissipated in the shunt resistances when the vortex moves 
with constant velocity. The result of this procedure for 
a square array is^ 

and for a triangular array (cf. Appendix A), 
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1] can also be obtained directly from the calculated /— {V) 
characteristics (cf. Appendix B). Our numerical results 
for different values of /3 in 8 x 8 arrays are shown in 
Table 2. This Table shows that r] varies approximately 
as /?^/^ at large (3 values, and differs considerably from 
the result (10). Indeed, at sufficiently large values of (3, 
the friction coefficient actually appears to be independent 
of the shunt resistance R. A similar result was also found 
previously in square arrays.''"^ 

A more accurate theoretical calculation of the fric- 
tional damping requires taking account of the loss of vor- 
tex energy to "spin-wave-like" excitations in the array. 
A continuum theory of this kind has been proposed by 
Geigenmiiller et aP . In Appendix C, we give a detailed, 
quasi-analytical theory for the friction coefficient, based 
on the loss of energy from a vortex moving with velocity 
V to spin wave modes. The model is more general than 
that of Ref. 9, in that it takes explicit account of the 
lattice structure of the array, so that a short-wavelength 
cutoff appears naturally, and allows for an arbitrary ex- 
ternal current source to excite the spin wave modes (for 
example, one arising from a high density of vortices). The 
final result for in a square array is 



and in a triangular array 



V 



iVo) 



tri 



87r2 ' 



(11) 



(12) 



where the dimensionless integral I' is given in Appendix 
C. 

It is sometimes useful to transform rj into an analogous 
expression for array resistivity. Using the force balance 

equation (8), and noting that the electric field has mag- 
nitude E = 27r(?i/(2e))n„w, where Uy = if / {^/Sa^) is the 
number of vortices per unit area, we obtain Ohm's Law 
in the form 

E = pJ, 

where the array resistivity is (after considerable simplifi- 
cation) 



32Tr^JR 
31' ■ 



(13) 



Table 3 shows the friction coefficient as calculated from 
the model of Appendix C. As can be seen, the resulting 
coefficient is strongly dependent upon /?, in agreement 



with our numerical "experiment." If. for example, we 
choose the scaled vortex velocity v = 1.0 (as defined in 
Eq. (32)), our analytic expression for r] is well approxi- 
mated in square array by the simple expression 



V « 0.89/3i/2(%).„ 



and in a triangular array by 



r] « 0.34/3^2(^0) 



(14) 



(15) 



The /3^/2 trend is consistent with the results of numerical 
calculations in Ref. 9 for square arrays and in the present 
work for triangular arrays, although the numerical coeffi- 
cient may differ by as much as a factor of two. The trend 
is also consistent with the experimental results of Ref. 8. 
Note that i) = 1 is a reasonable velocity to consider in 
this comparison, because larger velocities tend to trigger 
row-switching events in underdamped arrays. ^^^"^^^^^s 

Table 3 also lists the variation of ij with v at several 
values of /?. Evidently, r] is nearly independent of v at 
large v but goes to zero below a threshold value of or- 
der V = 0.2, where the damped pole in the integral (31) 
moves outside the first Brillouin zone of the triangular 
lattice. Once again, this agrees with the results found in 
the continuum theory of Ref. 9 for a square lattice. The 
results also agree quite well with our numerical results, 
as shown by a comparison of Tables 2 and 3(b). 

VI. ABSENCE OF BALLISTIC VORTEX MOTION 

A rather surprising result of our simulations is the 
absence of ballistic vortex motion. Such motion might 
have been expected, at least in the highly-underdamped 
regime; and there have been some experimental reports 
of such behavior^'^. In order to check for ballistic mo- 
tion, we apply a large bias current to the array in the 
flux flow regime, so that the vortex acquires a high ini- 
tial velocity, then we turn off the bias current. Since 
the effective mass of the vortex is presumably large in 
the high-/? regime, the vortex might be expected to move 
several lattice spacings because of its large initial kinetic 
energy, even after the driving current is removed. But 
from the calculated time-dependent voltage of the indi- 
vidual junctions, we find that the vortex travels at most 
through one primitive cell after the bias current is shut 
off, whatever its initial velocity, for all values of P con- 
sidered (0 < /3 < 1000). This is consistent with previous 
calculations''. 



6 



It is of some interest to compare our results with those 
of Ref. 10. This paper considers the triangular array 
in a continuum approximation and concludes that a nar- 
row vortex velocity window exists where ballistic motio- 
nis possible. It is suggested that this window extends 
from the vortex depinning current to roughly twice that 
current. We can envision two possible reasons why we do 
not see this ballistic regime in our own calculations. The 
first is that the depinning current is rather dependent 
on lattice size, typically being larger for the smaller lat- 
tices. For our size regime, this may narrow the window 
nearly to zero. In addition, the vortex velocity is not 
constant just above the depinning current, but instead 
is quite time-dependent, because of the periodic pinning 
potential. This time dependence is not considered in the 
model of ref. 10, which assumes a constant vortex ve- 
locity in estimating the width of the ballistic "window." 
Thus, the effects of this time-dependence could possibly 
also suppress this window. 

In view of these results, the explanation for the bal- 
listic motion which is observed in experiments seems un- 
clear. No numerical calculation has yet found such mo- 
tion from classical equations. Conceivably, the ballistic 
regime arises when quantum effects reduce the vortex 
friction coefficient below classical predictions, but this 
remains to be proven. 



ramctcrs where quantum corrections might need to be 
included in these calculations. In a naive picture, such 
corrections would start to matter when the characteris- 
tic charging energy (2e)^/C becomes comparable to the 
Josephson energy ?l7c/(2e). This condition gives 

This can be translated into a condition on the lattice 
resistivity, using Eqs. (12), (13), and (15), with the result 

^°^^((2^)' (1^) 

where the constant of proportionality is approximately 
1.1. This result suggests that quantum corrections might 
become important when the resistance per square of this 
two-dimensional network is comparable to the "quantum 
of resistance" /i/(2e)^ per imit frustration. Of course, 
this naive estimate does not take account of such obvious 
corrections as vortex- vortex interactions. It is amusing 
to note that several groups have reported evidence (both 
experimental and theoretical^*) for a superconductor-to- 
insulator transition in quasi-two-dimensional supercon- 
ductors in a magnetic field at a resistance per square of 
order /i/(2e)^; this transition is generally attributed to 
disorder effects, and thus may be unrelated to the sim- 
ple criterion for arrays mentioned above. Thus a detailed 
calculation of quantum effects on vortex motion remains 
an important problem for future study. 



VII. DISCUSSION AND CONCLUSIONS 

We have simulated the dynamical response of triangu- 
lar .Josephson junction arrays, using a classical model of 
resistively- and capacitively-shunted Josephson junctions 
described by coupled second-order differential equations. 
In the flux flow regime, the dynamical response of the 
network, including the time-dependent voltage, is well 
described in terms of vortex degrees of freedom. The vor- 
tices, however, experience higher viscous damping than 
predicted on the basis of a simplifled model, and in appar- 
ent contrast to experiment do not exhibit ballistic motion 
at any bias current we have investigated. The damping is 
reasonably well described, however, by a model which de- 
scribes loss of vortex energy to plasma (or "phase wave" ) 
oscillations in the Josephson network. 

It is of interest to make a crude estimate of the pa- 
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APPENDIX A: SIMPLE ESTIMATE OF VORTEX 
FRICTION COEFFICIENT 



In this Appendix, wc present a simple estimation of 
the friction cocfficiont for a single vortex moving in a 
triangular array. If we assume that such a vortex moves 
from the center of one plaquette to the center of a nearest- 
neighbor plaquette, it must cross one junction. Since the 
change in phase difference is 27r/3 when the vortex crosses 
the junction, the average voltage across the junction is 

2 



(y^ _ ^ 27r/3 _ h 27r/3 



h TT 

- ■ V, 



2e At 2e 2^a/v \/3 2e a ' 

where At is the time required for the vortex to move from 
one plaquette center to the next, a is the lattice constant, 
and V is the timc-avcragcd vortex velocity. 

Wc next compute the effective frictional coefficient by 
equating the resistively dissipated energy to that ex- 
pected for a particle moving in a viscous medium. Now 
the effective resistance between two nearest neighbor 
grains is defined as the voltage drop which is produced 
when a unit current is injected into one such grain and 
extracted from the other. Since there are six nearest 
neighbors in a triangular lattice, the eS'ective resistance 
in an infinite triangular array is R/3, where R is the 
single-junction resistance. Equating the resistively dis- 
sipated power to the frictional losses, we obtain 

1, ^ 2 3 4 ?i 2 ,^.2 2 

This implies that the vortex frictional coefficient in an 
infinite triangular array is 



(*)«=2.(A)».(f^±. 
A similar calculation for a square array gives^'^: 



(18) 



APPENDIX B: EXTRACTION OF r? FROM / - {V) 
CHARACTERISTICS 



In the flux flow region, the time- and space-averaged 
voltage {V) across an M x array is approximately pro- 
portional to the bias current /. We define a dimensionless 
proportionality coefficient 7 by 



{V) = -fNRI, 

where N is the number of junctions along the direction of 
the bias current. If we assume periodic transverse bound- 
ary conditions and [101] current injection, a complete cir- 
cuit of a vortex around the array produces a phase change 
of 27r across the array in the direction parallel to the cur- 
rent injection. The Josephson relation then implies that 

^ ' 2e^dt' 2eMa^' 
where a is the lattice constant and v is the transverse 
vortex velocity. Since the friction coefficient rjtri is re- 
lated to V by Eq. (8), we obtain, on combining the above 
relations. 



rjtri 



-yMN \2e 



1 

2R 



(20) 



APPENDIX C: ANALYTICAL MODEL FOR 
VORTEX FRICTION COEFFICIENT 



We consider a d-dimensional periodic network of 
RCSJ's, assuming zero shunt capacitance to ground, and 
also assuming that the shunt resistance, junction critical 
current, and shunt capacitance all vanish except between 
nearest neighbors, for which they take the values R, Ic, 
and C respectively. 

With these assumptions, the equations of motion for 
the phases may be written in the form 

h 

2e~ ^^'^ ' 2aR 



^'-C ^ hj + ^ 13 ^ij +Ic^ Sin(^ij) = li-ext, 



(21) 

where the sums run over the nearest neighbors to the 
grain i and (/>jj = (t)i — (jjj . 

Wc will calculate the losses produced by an externally 
applied current due to excitation of "spin waves", i. e. 
small-amplitude phase fluctuations. Thus, we expand the 
sine- function as sin(0j — (^j) w 0^ — and, within that 
assiimption, calculate the losses coming from an arbitrary 
Ii;ext{t)- With the introduction of the Fourier transforms 
(/)(k, w) and /ea:t(k, w), we can transform (21) into the 
form 
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/e,.t(k,w)A(k) 



T inuj 

-^c 2eR 2e 



(22) 



Here 



t(k) = 5^(l-exp(ik.R)), 



the sum runs over the set of nearest neighbor lattice vec- 
tors R, and the allowed k values run over the first Bril- 
louin zone of the grain lattice. 

The energy dissipation in the ij*^ bond in the time 
interval [-T, T] is: 



+T 

/^Eij = / I,^{t)Vij{t)dt, 

T 



or, in Fourier transform, 

/oo />oo 
dio I dio' hj{uj)V*j{J)A{uj,J), (23) 
-OO J —OO 

where 

A(.,.') = ^^*^. 

TT U — U' 

The total energy loss Ai^^ot = S<ij> . Using the 
Josephson relation between phase and voltage, and the 
equations of motion (21) in the small phase difference 

approximation, and making the relevant Fourier trans- 
forms, we can finally express the total energy loss as 



AS 



tot 



2eN 



x^^^frl^^W^-^C-.-'). (24) 

-"c "T 2eR 2e 

where TV is the number of grains in the lattice, and the 
sum runs over the first Brillouin zone. 

This result is valid for an arbitrary external current 
source. We now specialize to a vortex moving with ve- 
locity V. According to GeigenmiiUer et al, such a vortex 
traveling in the y direction has associated with it a charge 
density 

Q^(x, t) = C^v-^2Trd{x)6iy - vt). 

The corresponding space and time Fourier transform is: 

Q^(k,a;) = {27rf/^C^v{~ik,)S{u> - vky). (25) 

Of course, this charge density was derived for a vortex 
moving in a continuum superconductor, and cannot be 
exactly correct for a superconducting array. However, 



by imposing the additional requirement that Q^(k,CL') 
vanish for k outside the first Brillouin zone, we produce 
an approximate charge density which is properly discrete. 

The external current corresponding to this charge dis- 
tribution is: 

/ext(k,a;) = -iu}Q^{]s.,u}). 

We substitute this into Eq. (24) and convert the sum 
over k into an integral, with the help of a factor S'/(47r^), 
where S is the area of the array. The integrand vanishes 
unless uj = uj' . Next, wc explicitly evaluate the real part, 
using the fact that A{lu,uj) = ^, and dividing by 2T. The 
two integrals over frequency can be done immediately, 
since they involve delta-functions, and the final result for 
the time rate of energy loss from the vortex into the spin 
wave modes reduces to 



dE o 



(26) 



where rj is an effective vortex friction coefficient . The 
form (26) properly corresponds to a frictional force of 
the form —rjv, since the rate of energy loss is the dot 
product of the force with the velocity. 

The expression for r] can be written most compactly 
by using the dimensionless variables k'^ = aki, i = x 
or y, where a is the bond length. We also introduce a 
dimensionless vortex velocity 



v/{aoJo), 



(27) 



where luq = \/2eT^hC is the Josephson plasma fre- 
quency. After some algebra, we obtain Eqs. (11) and 
(12), respectively, for 77 in a square array and a trian- 
gular array. In both cases the dimensionless integral /' 
takes the form 



r'(«,/3)= \ 

Jb.z. 



dk^dJzy 



1 



(28) 



and (?7o)sg and {rio)tri are given by Eq. (9) and Eq. (10). 

The integral for both lattices runs over the scaled first 
Brillouin zone of the array (defined by taking the bond 
length a = 1). 
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FIGURE CAPTIONS 



the periodic boundary conditions are represented 
by repeating the N x N lattice. 



1. Schematic diagram of an 8 x 8 triangular Joseph- 
son junction array. Each intersection represents a 
superconducting grain, which is connected to its 
six nearest neighbors by Josephson coupling, (a) 
and (b) correspond respectively to [101] and [211] 
current injection direction. Free boundary condi- 
tions are used in the direction of current injection, 
while periodic boundary conditions are used in the 
transverse direction. 

2. Schematic illustration of a triangular plaquette of 
Josephson junctions at zero magnetic field, sub- 
jected to an injected current / as shown. The crit- 
ical current for this arrangement is 1.7602/c. 

3. / — {V) characteristics for 8 x 8 triangular arrays 
at / = 1/128 at two different (3 values with current 
injected in the [101] direction: (a) /3 = 0; (b) /? = 
10. The insets are enlargements of the flux flow 
regime. 

4. 7 — {V) characteristics for overdamped 8x8 tri- 
angular arrays {(3 = 0) at / = 1/128 with current 
injected in the [211] direction. The inset is the en- 
largement of the flux flow regime. 

5. Time dependent voltage traces and vortex motion 
path in the array for 8 x 8 overdamped arrays 
(/? = 0) at / = 1/128 with current injected in 
the [101] direction for two different applied cur- 
rents, to = Ti/(2eRIc) is a natural unit of time. 
The bias current and time-averaged voltages are 
(a) I/Ic = 0.2, {V)/NRIc = 2.39 x lO'^; (b) 
I/h = 1-0, {V)INRh = 1.91 X 10-2. (c) Path 
of vortex motion in this array; the disks represent 
successive positions of the vortex in the array. 

6. Time dependent voltage traces and vortex trajec- 
tories in an 8 X 8 overdamped array (/? = 0) at / = 
1/128, for three different values of current applied 
in the [211] direction, = h/{2eRIc) is the natural 
unit of time. The bias currents and time-averaged 
voltages are (a) I/h = 0.4, (V) /NRI^ = 4.298 x 
10-3; (b) I/Ic = 0.9, {V)/NRIc = 1.123 X 10"^; 
and (c) I/Ic = 1.3,{V)/NRIc = 1.978 x IQ-^. 
(d), (e), and (f) show the vortex trajectories cor- 
responding to (a), (b), and (c). We draw the tra- 
jectories in a "repeated lattice scheme" in which 
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TABLE CAPTIONS 



Table 3 (a) 



1 . Numerical values of the critical current for a,n NxN 
triangular array at / = 1/{2N^), for [101] current 
injection. 

2. Numerical values of 77 as a function of /? in an 8x8 
triangular array. 7 is estimated from the numerical 
I — {V) characteristics of the corresponding arrays. 
T] is calculated from Eq. (20) in Appendix B, and is 
estimated by drawing a straight line by eye through 
the / — (F) characteristic in the flux-flow regime. 

3. Numerical values of t] as obtained from the semi- 
analytical theory of Appendix C at several values 
of P and V for square and triangular arrays. The 
values are obtained by use of Eqs. (33) - (35), car- 
rying out the integral numerically. 



Table 1 



N 


8 


12 


16 


24 


Id/Ic 


0.090 


0.063 


0.054 


0.048 



Table 2 



p 


7 


v/no 





2.73 X 10"^ 


0.572 


1 


2.13 X 10-2 


0.732 


10 


8.57 X 10"* 


1.82 


50 


4.09 X 10"^ 


3.82 


100 


2.78 X 10"^ 


5.63 


225 


2.07 X 10"* 


7.54 


400 


1.47 X 10"* 


10.7 











1 


10 


50 


100 


225 


400 




0.2 


0.040 


0.061 


0.064 


0.064 


0.064 


0.064 




0.5 


0.67 


3.16 


7.67 


11.0 


16.7 


22.4 


V 


1.0 


1.14 


3.25 


6.75 


9.33 


13.7 


18.2 




2.0 


1.38 


2.82 


5.10 


6.83 


10.1 


14.1 



Table 3 (b) 







/3 




1 


10 


50 


100 


225 


400 




0.2 


0.014 


0.025 


0.027 


0.027 


0.028 


0.028 




0.5 


0.20 


0.94 


2.32 


3.36 


5.13 


6.91 


V 


1.0 


0.38 


1.19 


2.55 


3.56 


5.25 


6.84 




2.0 


0.48 


1.06 


1.98 


2.66 


3.86 


5.06 



Yu et al. Fig. 1 




Yu at al. Fig. 2 
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Yu et al. Fig. 6 (e) 




Yu et al. Fig. 6 (f) 




